program main
!
!                                              *******
! 	        Programa "Resolução da Equação de Schrodinger Dependente do Tempo
!              	 em Sistema Quase-Unidimensional com o Método de Split Operator"
!
!                                              *******
!		        		   ** Versão 0.8 **
!			               Estefferson Torres Carlos
!               	             Orientador: Gil de Aquino Farias
!                                      Co-Orientador: Andrey Chaves
!
! *** Changelog ***
! * v0.8
! -- modificado para trabalhar com partícula rotacionando em torno de um ponto em trajetória circular
!
! * v0.7
! -- adicionado potencial de múltiplos poços/barreiras quadradas finitas
! -- agora no arquivo "psi-0.dt" e outros arquivos é escrito um número de colunas exatamente igual ao número de estados
!
! * v0.6
! -- adicionado a evolução de uma onda em um potencial que varia com o tempo (pulso luminoso);
!
! * v0.5
! -- adicionado o cálculo dos N primeiros autoestados;
! -- adicionado o cálculo das projeções de uma função dada sobre cada um dos auto-estados;
! -- adicionado um timer para contabilizar o tempo de execução do programa;
!
! * v0.4
! -- adicionado o cálculo dos 5 primeiros autoestados (fundamental e quatro excitados);
!
! * v0.3
! -- adicionado poço de potencial quadrado;
! -- adicionado o tempo complexo;
! -- adicionado cálculo da primeiro energia do sistema;
! -- substituído o contador temporal "j" por "it";
! -- função do potencial V transformada em subrotina;
!
! * v0.2
! -- adiconado o uso de massa variável;
! -- rotina do método de Gram-Schmidt disponível, mas não usada;
!
! * v0.1
! -- versão inicial
!
! *** Descrição do Problema ***
! A equação de Schrodinger dependente do tempo tem a seguinte forma:
!
! ih * d(psi(x,t))/dt = H*psi(x,t)
!
! A função psi(x,t) descreve o pacote de onda em t=0 e possui a seguinte forma:
!
! psi(x,t) = exp(-x**2/gamma)
!
! Onde x é a espaço no eixo x e gamma é uma constante real arbitrária
!
! Depois da resolução anlítica do problema e da discretização do mesmo para
! a resolução numérica, tem-se um sistema de equações lineares da seguinte forma:
!
! A(i)*eta(i-1) + B(i)*eta(i) + C(i)*eta(i-1) = theta(i)
!
! Neste programa teremos adicionados além da massa variável, o tempo complexo, o cáclculo de autoestados
! e das energias destes autoestados.
! *** Sem Variáveis Implícitas
	implicit none

! *** DECLARAÇÃO DE VARIÁVEIS E CONSTANTES ***

! * CONSTANTES INTEIRAS
! 'nangle' é o número de subintervalos em que o sistema está dividido na direção x
! 'nt' é o número de iterações temporais
	integer, parameter :: nangle = 10001, NSTATES = 4, nt = 20000
	integer, parameter :: outunit=45, outunit2=46

! * VARIÁVEIS CARACTERE
! Variáveis usadas para gerar .dat para criar uma animação
! 'fn' recebe o nome dos arquivos .dat a serem usados na animação
! 'addleadzeros' é a função que adiciona algarismos '0' no número
! 'fn_number' é a numeração a ser adicionada ao nome do arquivo
	character(len=70) :: fn
	character(len=20) :: addleadzeros, fn_number

! * VARIÁVEIS INTEIRAS
! 'i' e a variável de incremento a ser usadas nos laços em geral
! 'j' é a variável de incremento a ser usada no laço mais externo de iteração temporal
!	- nota: para programas em 2D ou 3D basta substituir i por ix, iy e iz; e j por it
! 'imr' determina se o tempo será complexo (=1) ou real (=0)
! 'is' determina o número de estados a serem trabalhados
	integer i, it, imr, is

! * VARIÁVEIS REAIS
! 'hc' é a constante h cortado (dada em meV*fs)
! 'a0' é o raio de Bohr (dado em Angstrom)
! 'Ry' é a energia de Rydberg (dada em meV)
! 'am' é a massa do elétron (em elétron-massa)
! 'am1' e 'am2' são valores de massa expressos em função de 'am'
! 'ec' é a carga do elétron
! 'pi' é o número pi
! 'dt' é o intervalo de tempo (dado em fs - femtosegundos)
	double precision hc, a0, Ry, am, am1, am2, ec, pi

! 'radius' é o raio da circunferência na qual a partícula rotaciona
! 'langlet' é a largura de todo o sistema na direção x (em angstrom)
! 'lx' é a largura do poço na direção x (em angstrom)
! 'dangle' é o tamanho do passo na diração x (em angstrom)
! 'gamatheta' é a largura da gaussiana inicial do pacote de onda
! 'normal' é o fator de normalização da função de onda (somente da inicial?)
! 'x0' é a posição inicial do pacote de onda
! 'norm' é o fator de normalização
! 'mass' é a massa em unidades da massa do elétron
! 'F' é a amplitude do campo elétrico que será utilizado em alguns cálculos
! 'potwidth' é a lagura do poço de potencial quadrado
! 'Ef' é a energia do fóton que entra no sistema
	double precision radius, langlet, dangle, gamatheta, x0, xsum, xavg, psisum, F, potwidth, potwidth1, Ef
!	double precision elaptime, time1
	real etime          ! Declare the type of etime()
	real elapsedtime(2) ! For receiving user and system time
	real totaltime      ! For receiving total time

! 'angle(i)' é a coordenada espacial do eixo x
! 'V(i)' é o potencial
! 'mass(i)' é a massa
! 'E(i)' é a energia
	double precision, dimension(nangle) :: angle, V
	double precision, dimension(0:nangle+1) :: mass
	double precision, dimension(NSTATES) :: E!, proj
	double precision asum1, asum2, asum3
	double complex, dimension(NSTATES) :: proj
!	real, dimension(2) :: tarray
!	real :: resulttime

! * VARIÁVEIS COMPLEXAS
! 'img' é o número imaginário puro "i"!
! 'eps' é o epsilon = i*h*dt/2m
! 'dt' é o tempo complexo (dado em fs - femtasum2 = asum2 + cdabs(psi(i,2))**2*dthetaosegundos)
	double complex :: img, dt, eps, cerror, cerf

! 'AL', 'BL' e 'CL' são coeficientes do sistema
! 'AR', 'BR' e 'CR' são coeficientes do sistema
! 'theta' é a matriz-coluna de coeficientes da direita
! 'eta' é o array com as soluções eta(x,t+dt)
! 'xsi' é o array com a função xsi(x,t+dt)
! 'psi' é o array com a função psi(x,t)
	double complex, dimension(nangle) :: AL, BL, CL, AR, BR, CR, theta, eta, xsi, cexpv
	double complex, dimension(nangle,NSTATES) :: psi, psieig

	common / constants / dangle, dt, img, Ry, a0, am, am1, am2, ec, pi, radius

! *** FIM DECLARAÇÃO DE VARIÁVEIS E CONSTANTES ***

! *** ABERTURA DAS UNIDADES DE GRAVAÇÃO DE DADOS ***
! Abrindo as unidades para gravação de dados
	open(unit=10,file='psi-0.dat')
	open(unit=20,file='psi-dt.dat')
	open(unit=30,file='xavg.dat')
	open(unit=40,file='abc.dat')
! AVISO: O NÚMERO 45 ESTÁ RESERVADO PARA ESCREVER OS ARQUIVOS .DAT PARA A ANIMAÇÃO
	open(unit=50,file='potential.dat')
	open(unit=51,file='mass.dat')
	open(unit=52,file='energies.dat')
	open(unit=60,file='psi-0-proj.dat')
	open(unit=61,file='psi-dt-proj.dat')
	open(unit=62,file='xavg-proj.dat')
!	open(unit=63,file='potential-proj.dat')
	open(unit=64,file='projections.dat')
	open(unit=65,file='field-proj.dat')
	open(unit=66,file='potwidth-proj.dat')
	open(unit=67,file='x-time-psi-dt.dat')
	open(unit=68,file='photon-energy-proj.dat')	

!	open(unit=70,file='psi-0-proj-vt.dat')
!	open(unit=71,file='psi-dt-proj-vt.dat')
!	open(unit=72,file='xavg-proj-vt.dat')
!	open(unit=73,file='potential-proj-vt.dat')
!	open(unit=74,file='projections-vt.dat')

! *** ATRIBUÇÃO DE VALORES ***

! Chama a subrotina que com informações sobre o programa
call proginfo

! * Constantes Complexas
! Valor da unidade imaginária
img = CMPLX(0.D0,1.D0)

imr = 1 ! = 0, propagação em tempo real
        ! = 1, propagação em tempo imaginário

! * Constantes Reais
a0 = 0.5292D0
Ry = 13.6057D3
hc = 1.054589D3/1.602189D0
pi = 4.D0*datan(1.D0)
ec = 1.60217657D-19
! Massas efetivas
am = 0.067D0
!am = 1.D0
am1 = 1.D0*am
am2 = 1.D0*am

! Constantes temporais
dt = 0.1D0
if(imr.eq.1) dt = -dt*img
!write(*,*) dt

! Valor do raio da circunferência
radius = 100.D0/a0

! Valor das larguras do sistema e do passo dangle
langlet = 2*pi
dangle = langlet/(nangle-1)

!write(*,225) langlet, dangle

! Valor de angle(i)
x0 = 0.D0/a0

do i=1,nangle
!	angle(i) = (-langlet/2.D0) + (i-1)*dangle
        angle(i) = (i-1)*dangle
!	write(*,190) i, angle(i)
!	pause
enddo

!stop

! * PACOTE DE ONDA INICIAL
! Valor de gama
gamatheta = 500.D0/a0

! Chamando a subrotina initpsi para escrever o pacote de onde inicial
call initpsi(x0,gamatheta,angle,psi,NSTATES,imr)

! Escrevendo o pacote de onda inicial no arquivo psi-0.dat
do i=1,nangle
	write(10,220) angle(i), angle(i)/pi, ((cdabs(psi(i,is))**2),is=1,NSTATES)
!	write(*,220) angle(i), ((cdabs(psi(i,is))**2),is=1,NSTATES)
enddo

! * POTENCIAL
! Valor de Potencial para todo o programa
potwidth1 = 400.D0/a0

! Chama potcalc para calcular o potencial
potwidth = potwidth1
!write(*,*) potwidth*a0

! Valor do Campo Elétrico F (em KV/cm)
F = 20.D0

!write(*,*) radius, radius*a0
!stop

call potcalc(angle,V,nangle,0,F,radius,0.D0)

! Escreve V(i) no arquivo potencial.dat
do i=1,nangle
	write(50,201) angle(i)/pi, V(i)*Ry
!	write(*,216) angle(i), angle(i)/pi, V(i)*Ry
enddo

!stop

! * MASSA
! Chama masscalc para calcular mass(i)
call masscalc(angle,mass)

! Escreve mass(i) no arquivo mass.dat
do i=1,nangle
	write(51,216) angle(i), angle(i)/pi, mass(i)
!	write(*,216) angle(i), angle(i)/pi, mass(i)
enddo

!stop

! * Variáveis Complexas
! Valor de eps
	!eps = (img*dt)/(hc)
	eps = -(img*dt*Ry)/(hc)

! *** FIM ATRIBUÇÃO DE VALORES ***

! * COEFICIENTES DAS MATRIZES TRIDIAGONAIS
! Laço para calcular os valores de A(i), B(i) e C(i)
! !!!ATENÇÃO!!!: O problema não com os coeficientes, pois os valores aqui e no programa do Andrey são praticamente
! os mesmos! O erro está em outro lugar!!!
! O erro deve estar em algum dos trechos seguintes, ou então está relacionado à constantes que não estão colocadas
! com o valor correto.
! O erro deve estar em algum termo que deveria ter sido alterado devido à inclusão do tempo imaginário.
do i=1,nangle
	cexpv(i) = cdexp(eps*V(i)/2.D0)
! Coeficientes do lado esquerdo
	AL(i) = - (eps * (-1.D0/dangle**2/4.D0) * ((1.D0/mass(i-1)) + (1.D0/mass(i))))
	BL(i) = (1.D0 - (eps * (1.D0/dangle**2/4.D0) * ((1.D0/mass(i-1)) + (2.D0/mass(i)) + (1.D0/mass(i+1)))))
	CL(i) = - (eps * (-1.D0/dangle**2/4.D0) * ((1.D0/mass(i)) + (1.D0/mass(i+1))))
! Coeficientes do lado direito
	AR(i) = (eps * (-1.D0/dangle**2/4.D0) * ((1.D0/mass(i-1)) + (1.D0/mass(i))))
	BR(i) = (1.D0 + (eps * (1.D0/dangle**2/4.D0) * ((1.D0/mass(i-1)) + (2.D0/mass(i)) + (1.D0/mass(i+1)))))
	CR(i) = (eps * (-1.D0/dangle**2/4.D0) * ((1.D0/mass(i)) + (1.D0/mass(i+1))))
enddo

do i=1,nangle
	write(40,225) cdabs(AL(i)), (BL(i)), cdabs(CL(i)), cdabs(AR(i)), cdabs(BR(i)), cdabs(CR(i))
!	write(*,225) cdabs(AL(i)), cdabs(BL(i)), cdabs(CL(i)), cdabs(AR(i)), cdabs(BR(i)), cdabs(CR(i))
!	write(*,225) AL(i), BL(i), CL(i), AR(i), BR(i), CR(i)
enddo

!stop

! =========================================================================================
!					EVOLUÇÃO DA GAUSSIANA EM TEMPO REAL/IMAGINÁRIO
! =========================================================================================
write(*,*) "Evolução para encontrar os autoestados começou..."

! Começa aqui o laço para calcular o valor de psi(i,is)
do it=1,nt
	do is=1,NSTATES
		do i=1,nangle
			xsi(i) = cexpv(i) * psi(i,is) ! Valor de xsi(i)
		enddo

! Valor de theta(i) - a coluna da direita
		theta(1) = BR(1)*xsi(1) + CR(1)*xsi(2)
		theta(nangle) = AR(nangle)*xsi(nangle-1) + BR(nangle)*xsi(nangle)

		do i=2,nangle-1
			theta(i) = AR(i)*xsi(i-1) + BR(i)*xsi(i) + CR(i)*xsi(i+1)
		enddo

! O valor de eta(i) será retornado pela subrotina tridag, por isso não será atribuido
! A subrotina TRIDAGX(a,b,c,r,u) é chamada para resolver o sistema
		call TRIDAGX(AL,BL,CL,theta,eta)

! Multiplica-se a função eta(i) por
		do i=1,nangle
			psi(i,is) = cexpv(i) * eta(i)
		enddo

! Normalização de psi a cada iteração temporal
		if(imr.eq.1)then
			psisum = 0.D0
			do i = 1,nangle
				psisum = psisum + cdabs(psi(i,is))**2*dangle
			enddo

			do i = 1,nangle
				psi(i,is) = psi(i,is)/dsqrt(psisum) !dividing by the normalizing factor
			enddo

			call GRAMSCHMIDT(psi,NSTATES)

			if(mod(it,10).eq.0.D0)then
				call ENERGY(mass,V,psi,E)
				if(isnan(E(is)))then
					E = 0.d0
                                        write(*,*) "A energia estourou! Saindo do programa..."
                                        stop
				endif
				write(52,181) it*cdabs(dt), E
				write(*,181) it*cdabs(dt), E
			endif
		endif

! Escrevendo a posição média da frente de onda
		if(imr.eq.0)then
			xsum = 0.D0
			do i = 1,nangle
				xsum = xsum + cdabs(psi(i,is))**2*angle(i)*dangle
			enddo
			xavg = xsum
			write(30,225) cdabs(it*dt),xavg*a0,xavg*a0*ec
		endif
	enddo ! fim do loop dos estados
enddo ! fim do loop temporal

! Escrevendo o pacote de onda ao final de todas as iterações temporais
do i=1,nangle
	write(20,220) angle(i), ((cdabs(psi(i,is))**2),is=1,NSTATES)!cdabs(psi(i,1))**2, cdabs(psi(i,2))**2, cdabs(psi(i,3))**2, cdabs(psi(i,4))**2, cdabs(psi(i,5))**2, cdabs(psi(i,6))**2
enddo

do is=1,NSTATES
	do i=1,nangle
		psieig(i,is) = psi(i,is)
	enddo
enddo

stop

write(*,*) "Evolução para encontrar os autoestados terminou..."

stop

! =========================================================================================
! 					EVOLUÇÃO EM TEMPO REAL DO PRIMEIRO AUTESTADO
! =========================================================================================
! Armazenando o valor de psi(i,is) na variável psieig(i,is). A função psieig servirá para
! armazenar os valores das funções dos autoestados para que depois seja feita a projeção
! de psi sobre cada um dos autoestados.

write(*,*) "Evolução em tempo real do primeiro autoestado começou..."

! Agora, a função psi(i,1) = psieig(i,1) será a função inicial, e será evoluída em tempo real,
! re-usando os mesmos algorítimos usados na primeira parte do programa.
! Para garantir que psi(i,1) será igual a psieig(i,1):
do i=1,nangle
	psi(i,1) = psieig(i,1)
	write(60,220) angle(i), cdabs(psi(i,1))**2
enddo

! * Atribuições de Variáveis
imr = 0 !para que o sistema evolua em tempo real
dt = 0.1D0 !para que o tempo seja real

! Valor de eps
eps = -(img*dt*Ry)/(hc)

! * COEFICIENTES DAS MATRIZES TRIDIAGONAIS
! Laço para calcular os valores de A(i), B(i) e C(i)
do i=1,nangle
! Coeficientes do lado esquerdo
	AL(i) = - (eps * (-1.D0/dangle**2/4.D0) * ((1.D0/mass(i-1)) + (1.D0/mass(i))))
	BL(i) = (1.D0 - (eps * (1.D0/dangle**2/4.D0) * ((1.D0/mass(i-1)) + (2.D0/mass(i)) + (1.D0/mass(i+1)))))
	CL(i) = - (eps * (-1.D0/dangle**2/4.D0) * ((1.D0/mass(i)) + (1.D0/mass(i+1))))
! Coeficientes do lado direito
	AR(i) = (eps * (-1.D0/dangle**2/4.D0) * ((1.D0/mass(i-1)) + (1.D0/mass(i))))
	BR(i) = (1.D0 + (eps * (1.D0/dangle**2/4.D0) * ((1.D0/mass(i-1)) + (2.D0/mass(i)) + (1.D0/mass(i+1)))))
	CR(i) = (eps * (-1.D0/dangle**2/4.D0) * ((1.D0/mass(i)) + (1.D0/mass(i+1))))
enddo

potwidth = potwidth1
! Valor da energia do fóton Ef (em meV)
Ef = 4.D0
!write(*,*) potwidth*a0
!pause

! * Evoluindo o novo pacote de onda inicial psi(i,1) no tempo real
do it=1,nt
! * POTENCIAL
! Chama potcalc para calcular o potencial
		call potcalc(angle,V,nangle,it,F,potwidth,Ef)
! Escreve V(i) no arquivo potencial.dat (sei que gasta tempo, mas vou deixar assim)
!		do i=1,nangle
!			write(63,200) angle(i), V(i)*Ry
!		enddo
! Escrevendo o termo exponencial com o potencial
		do i=1,nangle
			cexpv(i) = cdexp(eps*V(i)/2.D0)
		enddo

		do i=1,nangle
			xsi(i) = cexpv(i) * psi(i,1) ! Valor de xsi(i)
		enddo

! Valor de theta(i) - a coluna da direita
		theta(1) = BR(1)*xsi(1) + CR(1)*xsi(2)
		theta(nangle) = AR(nangle)*xsi(nangle-1) + BR(nangle)*xsi(nangle)

		do i=2,nangle-1
			theta(i) = AR(i)*xsi(i-1) + BR(i)*xsi(i) + CR(i)*xsi(i+1)
		enddo

! O valor de eta(i) será retornado pela subrotina tridag, por isso não será atribuido
! A subrotina TRIDAGX(a,b,c,r,u) é chamada para resolver o sistema
		call TRIDAGX(AL,BL,CL,theta,eta)

! Multiplica-se a função eta(i) poasum2 = asum2 + cdabs(psi(i,2))**2*dthetar
		do i=1,nangle
			psi(i,1) = cexpv(i) * eta(i)
		enddo

! Normalização de psi a cada iteração temporal
! NOTA: 24/10/2012 - Retirei a normalização para ver se era isto que fazia com que
! a soma das projeções fosse maior que 1, mas o resultado é praticamente o mesmo.
		if(imr.eq.0)then
			psisum = 0.D0
			do i = 1,nangle
				psisum = psisum + cdabs(psi(i,1))**2*dangle
			enddo

			do i = 1,nangle
				psi(i,1) = psi(i,1)/dsqrt(psisum) !dividing by the normalizing factor
			enddo
		endif

		if(imr.eq.1)then
			call GRAMSCHMIDT(psi,NSTATES)

			if(mod(it,10).eq.0.D0)then
				call ENERGY(mass,V,psi,E)
				write(52,180) it*cdabs(dt), E
				write(*,180) it*cdabs(dt), E
			endif
		endif

! Escrevendo a posição média da frente de onda
		if(imr.eq.0)then
			xsum = 0.D0
			do i = 1,nangle
				xsum = xsum + cdabs(psi(i,1))**2*angle(i)*dangle
			enddo
			xavg = xsum
			write(62,225) cdabs(it*dt),xavg*a0,xavg*a0*ec
		endif

! Chamando a subrotina projcalc para calcular as projeções de psi sobre psieig
		if(mod(it,10).eq.0.D0)then
			call projcalc(proj,psi,psieig,angle,nangle,NSTATES)
! Escrevendo o valor do projeção no arquivo 'projections.dat'
			write(64,220) it*cdabs(dt), cdabs(proj)**2
            if(mod(it,1000 ).eq.0.D0)then
    			write(*,220) it*cdabs(dt), cdabs(proj)**2
    		endif
		endif

! If para escrever um video/pot-000000000000x.dat, que será usado para criar uma animação
! com a evolução do pacote de onda. Por enquanto apenas para o tempo real.
        if(imr.eq.0)then
            if(mod(it,50).eq.0.D0)then
                ! create the directory if it doesn't exist
                if(it.eq.1)then
                        call SYSTEM('if [ ! -d "potvideo" ]; then mkdir "potvideo"; fi; cd "potvideo"; rm *.dat 2> /dev/null')
                endif
                ! build filename -- i.dat
                fn_number = addleadzeros(it,20)
                write(fn,fmt='(a)') 'potvideo/' // 'pot-' // fn_number // '.dat'
                ! open it with a fixed unit number
                open(unit=outunit,file=fn, form='formatted')
                ! write something
!               write(*,*) 'Escrevendo o potencial em', fn
                do i = 1,nangle
                    write(outunit,220) (angle(i))/10.D0, V(i)*Ry
                enddo
                ! close it
                close(outunit)
            endif
        endif

! If para escrever um video/psi-000000000000x.dat, que será usado para criar uma animação
! com a evolução do pacote de onda. Por enquanto apenas para o tempo real.
        if(imr.eq.0)then
            if(mod(it,50).eq.0.D0)then
                ! create the directory if it doesn't exist
                if(it.eq.1)then
                        call SYSTEM('if [ ! -d "psivideo" ]; then mkdir "psivideo"; fi; cd "psivideo"; rm *.dat 2> /dev/null')
                endif
                ! build filename -- i.dat
                fn_number = addleadzeros(it,20)
                write(fn,fmt='(a)') 'psivideo/' // 'psi-' // fn_number // '.dat'
                ! open it with a fixed unit number
                open(unit=outunit2,file=fn, form='formatted')
                ! write something
!               write(*,*) 'Escrevendo o potencial em', fn
                do i = 1,nangle
                    write(outunit2,220) (angle(i))/10.D0, cdabs(psi(i,1))**2
                enddo
                ! close it
                close(outunit)
            endif
        endif

! Escrevendo os valores de posição, tempo e função de onde evoluída para
! usar no mapa de cores, mostrando como a função varia no tempo
		if(mod(it,100).eq.0.D0)then
!			do i=1,nangle
			do i=3000,7000
				write(67,222) (angle(i))/10.D0, it*cdabs(dt), cdabs(psi(i,1))**2
			enddo
		endif

enddo ! fim do loop temporal

write(65,220) F, cdabs(proj)**2
write(66,220) potwidth*a0, cdabs(proj)**2
write(68,220) Ef, cdabs(proj)**2

do i=1,nangle
	write(61,220) angle(i), cdabs(psi(i,1))**2
enddo

write(*,*) "Evolução em tempo real do primeiro autoestado terminou..."

! ==========================================================================================
! 																					FORMATOS
! ==========================================================================================

100 FORMAT(I5,F16.3,20E16.6)
150 FORMAT(I5,2F16.3,20E16.6)
175 FORMAT(I5,F16.6)
180 FORMAT(10F16.6)
181 FORMAT(F16.6,50E16.6)
190 FORMAT(I5,6F16.6)
200 FORMAT(F16.3,2E16.6)
201 FORMAT(F16.6,2E16.6)
215 FORMAT(F16.6,E16.6)
216 FORMAT(2F16.9,E16.6)
220 FORMAT(F16.6,20E16.6E3)
222 FORMAT(2F16.6,1E16.6)
225 FORMAT(20E16.6)
250 FORMAT(2F16.3)
255 FORMAT(6F16.3)
300 FORMAT(I5,20E16.6)
400 FORMAT(2F16.3,2E16.6)

! Tempo para a execução do programa
	totaltime = etime(elapsedtime)
	print *, 'End: total=', totaltime, 'seg     user=', elapsedtime(1), 'seg    system=', elapsedtime(2)

!	call ETIME(tarray,resulttime)
!	write(*,*) "Elapsed time:", tarray(1), "seconds"
!	write(*,*) "system time:", tarray(2)
!	write(*,*) "total  time:", resulttime

end program main
